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A replicator equation with mutation processes is numerically studied. Without any mutations, 
two characteristics of the replicator dynamics are known: an exponential divergence of the 
dominance period, and hierarchical orderings of the attractors. A mutation introduces some 
new aspects: the emergence of structurally stable attractors, and chaotic itinerant behavior. In 
addition, it is reported that a neutral attractor can exist in the ji —* +0 region. 
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The replicator equation was initially proposed by May- 
nard Smith, EP and was developed thereafter to describe 
evolutionary dynamics (see e.g. HoibiauercP ) . ft is equiv- 
alent to the Lotka-Volterra equationjj-' and is now widely 
accepted as a basic model equation with applications 
ranging from ecosystems to other hierarchical network 
systems. 

Recently, unexpected rich behaviors have been found 
in the equations uW Of these, the emergence of an expo- 
nential time scale and of a complex attractor hierarchy 
are worth noting. The phenomenon of long-time domi- 
nance by a unique variable (species) and a chaotic tran- 
sition from one dominant species to another are generic 
features of the equation. The lifetime of the dominant 
species is found to diverge exponentially. This specific 
behavior is due to the heteroclinic cycles embedded in 
the replicator equation. Depending on one parameter, 
heteroclinic cycles can be hierarchically organized. 

On the other hand, the heteroclinic cycle has been 
considered to be unrealistic in the light of biological sys- 
tems. For example, a population size that decreases to 
the order of O(e -10 °) is considered to be unnatural in 
a real ecosystem. One remedy for this is to set a lower 
bound to the population size, as in the work of Tokita 
and Yasutomi.tiP That is, a species whose population 
drops below the given threshold must be removed from 
the system. As a result, the system in the end attains 
a stable distribution of species. However, as a conse- 
quence, the system loses its rich temporal behavior and 
many degrees of freedom. 

In the present letter, we propose another remedy: 
namely, to recover structural stability by introducing dif- 
fusion terms (see alsoO'lP). The diffusion terms can be 
identified as immigrations and mutations in an ecologi- 
cal system. We report on the emergence of heteroclinic 
chaos, chaotic itinerant phenomena, and neutral attrac- 
tors in the replicator equation with mutation. 
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Most studies on relatively large replicaffu; ^pmations 



deal with a random interaction matrix .LTtrB'LrEy This 
often makes it difficult to analyze the mechanism con- 
trolling the generic dynamical behavior. Here we intro- 
duce a symmetrical interaction matrix, and use a Lotka- 
Volterra-type equation. These two types of equations are 
mutually transformable. We use a Lotka-Volterra-type 
equation with seven degrees of freedom as the simplest 
model having more than one heteroclinic cycle. 

Xi = Xi [1 - Xi - a(xi+i + x i+2 + x i+i ) 

- b(xi-i + Xi-2 + Xi-i)] (f ) 

where Xi > 0(i = 0, 1, • • • , 6). There are two characteris- 
tic parameters a and 6, which satisfy the inequalities: 



a > 1 > b, 



-l>l-b 



(2) 



This equation is a natural extension of May's sys- 
tem.© We will discuss a family of replicator-type equa- 
tions which have several properties in common: f ) The 
R+(= {x\xi > 0(i = (),•••, 6)}) plane is kept in- 
variant. 2) There are saddle-type fixed points (e; = 
(0, 0, • • • , lj, • • • , 0)) at the boundaries and one inte- 
rior repeller-type fixed point, which is given by q{= 
i + 3a + 3b (l> !)■■•> !))■ 3) A null point gives a trivial fixed 
point (o= (0,0, •••,())). 

This equation has a symmetrical property in the sense 
that every saddle point has an equal number of incom- 
ing and outgoing directions of dimensionality one. For 
any z-th point, there are heteroclinic orbits from to 



= 1,2,4). 



As a result, there are seven saddle points and 21 hete- 
roclinic orbits in this system. Each incoming direction is 
correlated to one of the outgoing directions of the other 
saddle point (see Fig.l). All the initial points, except 
xq = ■ ■ ■ xq , will asymptotically converge to a network 
composed of the 21 heteroclinic orbits. _ 

As was first explicitly pointed out by Chawanyap 1 we 
also numerically find exponential divergence of the dom- 
inance period in the neighborhood of saddle points ej 



2 
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Fig. 1. A diagram of all possible heteroclinic orbits with seven 
saddle points symbolized as e^{k = 0,1, ■■•,6). Every saddle 
point has three incoming and three outgoing connections to/from 
other saddle points. For example, a saddle point ei has three 
outgoing connections to e2,ez and es and three incoming con- 
nections from erj,e4 and ee,. It is worth noting that e2, es and 
es compose another heteroclinic cycle of period 3, in addition 
to eo, e4 and eg. Because of this symmetry, each saddle point 
is taken to be equivalent. Thus, there exist three heteroclinic 
cycles which itinerate all the saddle points, but in a different 
order. They are drawn in different line styles. 



even with this symmetric equation (see Fig. 2, 4). Two 
different kinds of transition order from to ej are ob- 
tained. The transition order from ej to ej becomes cither 
chaotic (a), or periodic (b) with respect to initial state. 
Almost every case belongs to the chaotic case (a) . In the 
following, wc introduce a diffusion term into the equa- 
tion. By doing this, any heteroclinic cycle is removed 
from the system. Instead, we have the ruins of such cy- 
cles and the hoping dynamics among them. 



Fig. 2. Time vs. logiri, simulated with a = 1.2 and b = 0.9. The 
dominance periods in the neighborhood of each are gradually 
extended. When Xi is sufficiently close to the maximum value 
(i.e. unity), the growth dynamics of neighboring species are well 
approximated by the linear curves. That is, they are expressed 
by 1 - b and 1 - a for ^ log x i+j (j = 1, 2, 4) and ^ log x t -j (j = 
1,2,4) , respectively. 

Fig. 3. t - logx;: simulation run with the parameters a = 1.2, b = 
0.9,/i = 10~ 6 . It can be clearly observed that a lower bound 
exists to each population size, which is given approximately by 
-JU ~e- 12 - 2 . 

a — 1 

A possible diffusion process for mutation is introduced 
in the original replicator model as follows: 

Xi =Xi[l - Xi - a(x i+ i + x i+2 + x i+4 ) 

- b(xi-i + Xi-2 + Xi-i)] - 6fixi + /j,y^xj. (3) 

Here wc assume that there exists mutation from Xi to 
Xj in the ratio [i for all i and j. We call the diffusion 



Fig. 4. The n-th recurrence time for visiting a saddle point e; 
with parameters a = 1.2 and b = 0.9. 



process as mutation since it describes the flow from one 
population to the others, if we assume that the variable 
Xi as the population size of the species i and all species 
has the mutually transitionable genotypes. 

The mutation process naturally gives a lower bound 
to each population size (Fig. 3). Therefore, 
every equilibrium point (e^) on the peripherals van- 
ishes simultaneously. Consequently, heteroclinic cycles 
no longer exist. The exponential divergence of the dom- 
inance period is also suppressed. 

It is also worth noting that introducing mutation rates 
makes the system structually stable. Therefore, the nu- 
merical results are rather insensitive to the numerical 
method we used. At the same time, it is more reasonable 
to use a structually stable model to describe natural phe- 
nomena. Those two issues are great advantage to study 
the replicator equation with diffusion over other models. 

The equation is simulated numerically with the 4th 
order Runge-kutta method. Even the system is made 
structually stable, we solve an equation of a logarithm 
of each variable Xj. Doing this, we can study the very 
lower mutation rate regions. 

When the mutation rate is higher than t^t+IsT36)' 
the internal equilibrium point q becomes stable, which 
gives a unique fixed point in this system. Below the 
critical value, the q is destabilized and three limit cycles 
appear. Dominant species appear cyclically in each limit 
cycle with a fixed but different order(+l,+2 and +4). 
Each limit cycle is inversely characterized by this order. 
All the initial states in the phase space (R+) will be 
attracted to one of these limit cycles in this parameter 
region. 

By further decreasing the mutation rate, we see se- 
quential period bifurcation of each limit cycle, each hav- 
ing three quasi-periodic attractors(T 2 ) at some point. 
Each attractor holds the "dominant species recurrence 
order" which characterizes the original limit cycles found 
in the higher mutation-rate regime. It is difficult to ob- 
serve the higher order tori (T n ), and we instead observe 
the breakup of three T 2 , and the emergence of a strange 
attractor. The Lyapunov exponents are computed as a 
function of the mutation rate in order to quantify the 
observed route to chaos in Fig. 5. 

It is worth noting from Fig. 5 that there are two qual- 
itatively different chaotic attractors in the higher and 
lower mutation regimes: one with two positive Lyapunov 
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Fig. 5. log/i - A^: Lyapunov exponents computed for the param- 
eters a = 1.2 and b = 0.9 by varying the mutation rates. 



exponents, and one with a single positive exponent. 

We can interpret both structures of chaotic attractors 
as a combination of three local attractors connected by 
chaotic dynamics. These local attractors correspond to 
the heteroclinic cycles embedded in the original replica- 
tor system without mutations. 

The term chaotic itinerancy is used to describe a situ- 
ation where there is a natural separation of local attrac- 
tors and chaotic dynamics which connects those attrac- 
tors via a higher dimension subspace.E£P That is, a state 
itinerates chaotically among local attractors. 

If the local attractor itself is chaotic, it is difficult to 
distinguish whether the state is in or out of local attrac- 
tors. This case is typified by a chaotic attractor with 
two positive exponents. On the other hand, if the lo- 
cal attractor is not chaotic, it is possible to specify the 
state simply from the time-evolution of the size of the 
population. When a state is outside the local attractors, 
the population size changes rather chaotically. Indeed, 
we observe that the three local strange-attractors de- 
generate to quasi-periodic states in the lower mutation 
regime, where the dynamics connecting the local attrac- 
tors remains strange. Therefore, chaotic behavior and 
quasi-periodic behavior appear alternatively. 

This picture is well demonstrated by computing the 
local Lyapunov exponents. While the system stays in 
the local quasi-periodic state, the three largest expo- 
nents are computed as (0, 0, — ). While in the transition 
states, the exponents are computed as (+, 0, — ). Further, 
we estimate the local Kolmogorov-Sinai(KS) entropyEf 
to quantify the effective degrees of freedom at each mo- 
ment. The local KS entropy is estimated by the sum of 
the local Lyapunov exponents with positive values. The 
local Lyapunov exponents are computed from the local 
Jacobian of the dynamics. 

9 

(local)KSE = A„ (X q > 0, A g+ i < 0) (4) 
n=l 

The time series of the local KS entropy is plotted in 
Fig. 6, where intermittent bursts of the KS entropy are 
clearly observed. It is seen from this figure that the 
switching between local attractors is associated with a 
burst of the local KS entropy. Since the behavior of the 
local KS entropy is almost the same as that of the num- 
ber of the positive local exponents, we argue that this 



3 

switching behavior is followed by an increment in the 
number of effective degrees of freedom. From this char- 
acterization, we insist that the switching behavior found 
in this system should be named chaotic itinerancy. 



Fig. 6. t - localK S E: Local KS entropy is computed under the 
condition a = 1.2, b = 0.9 and fi = 10~ 6 . Where the local KS 
entropy is near zero, the orbit is in the neighborhood of one of 
the three local attractors. On the other hand, where the local 
KS entropy is large and positive, the orbit is in a transition state 
between local attractors. 



In the limit of fj, — » +0, the original heteroclinic be- 
havior is restored. Without any mutation terms, the 
system can show exponential growth of a single-species 
dominance period at u. = 0. On the other hand, we 
know that even a small mutation rate can remove the 
heteroclinic cycles (i.e. there emerges a lower boundary 
to a population size denoted by L^). Therefore, we ex- 
pect that fi = is a singular point. When we approach 
the point from the above (i.e. fx — > +0), the largest 
exponent will approach zero. This implies that the ex- 
pected dynamics becomes no more chaotic at this limit. 
Indeed, we see three different periodic behaviors which 
correspond to three heteroclinic cycles of the original 
replicator equation. Those periodic behaviors are thus 
indexed by a recurrent order of the dominant species(i.e. 
+ l,+2 and +4). It is worth noting that the motions 
with the same recurrence order densly construct a par- 
tially disconnected torus. That is, the limiting behavior 
is constrained on the torus. In the following we analyze 
the situation. 

Assuming that the mutation rate is inifinitesimally 
small, we can approximate the dynamics by some lin- 
earized equations. The limiting behavior means that the 
orbit is well charaterized by the transition between sad- 
dle points, and the approximation holds good enough 
around saddle points. There, the growth term A loga^ is 
well approximated by 1 — a, 1 — b and for the incoming, 
outgoing, and the dominant species, respectively. It is 
numerically found that each species becomes dominant 
twice, and diminishes (to the lowest order) four times 
during one period of time (see Fig. 7). This is a necessary 
feature in order to keep the periodicity of this linearized 
equation. It can be rewritten more formally in terms of 
successive dominance periods of given species, which are 
denoted by t\ and (see Fig. 7 for the meaning). 

(1 - 6)ta + (1 - a)tx < 0, (5) 
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(1 - b)h + (1 - a)t 2 < 0, (6) 

h + t * = -jHy ™ 



Fig. 7. With sufficiently small /i, the growth factor log 2;^ is 
approximated by 1 — a, 1 — 6 and for the incoming, outgoing, 
and the dominant species, respectively. 



The first and second inquality is derived by the fact 
that each species has to be at the lower boundary before 
it starts to dominate the population again. The third 
term is derived in order to hold the periodicity. That is, 
the sum of t\ and t 2 should be conserved by the dynam- 
ics, but its ratio |j can be redundant under the following 
condition: 



Within this range, absolute values of t\ and t 2 can be 
freely determined. This degree of freedom gives a neutral 
direction of the periodic state. That is, the periodic state 
can exist infinitely many. 

The above inequalities further set limits on where 
those periodic states can be found. We already know 
that there are three kinds of periodic cycle with respect 
to the recurrence order of the dominat species (+1,-1-2 
and +4). Each of those three cycles can exist infinitely 
and being found on the unique torus. However, they 
have to coexist on the torus and the boundaries between 
different cycles are given by the inequalities. Namely, 
(a — 1)/(1 — b) gives the maximal axial length of the 
partial torus size, where each cycle can occupy. There- 
fore, the torus is fragmented into three parts, each cor- 
responds to different cycle state ( in Fig. 8 (a)). 

Even with a small mutation rate, we see that the dy- 
namics takes off from the torus. The neutral directions 
discussed above does not exist anymore. The dynamics 
with the finite mutation rate is illuastrated in Fig.8(b). 
The periodic orbit gradually shifts to the edges of the re- 
gions, and finally switches off to other regions via chaotic 
motions. The original periodic behaviors now become 
quasi-periodic. 

In this letter, we have shown how mutation dynamics 
changes the behavior of the original replicator system. 
In particular, a chaotic itinerancy is noted. We under- 
stand phenomenologically that the dynamic behavior of 
this system is well described by the composition of local 
attractors, and hoping between them, where the local 
attractors are ruins of heteroclinic cycles. 

This system is made symmetric so that each saddle 



(a) (b) 

Fig. 8. A schematic drawing of the dynamics at fi — » +0 (a) and 
at fi = +e (b). An orbit asymptotically converges to a cycle 
on a transverse direction of this torus at fi — > +0. But with a 
small mutation rate, the orbit will no longer stay in one portion 
of the torus, but switches from one to the other. The switching 
behavior appears to be chaotic. 



point has an equal number of incoming and outgoing 
connections to other saddle points. Understanding a bi- 
furcation diagram and the limiting behavior was possible 
due to the symmetry of the system. 

What we have observed in this system is not restricted 
to the present symmetrical cases. Replicator systems 
with partially disordered bimatrix have been studied, 
and the corresponding behavior has been reported (tP). 
A direct extension from this system is to study a series of 
symmetric equations with more than three local attrac- 
tors. We would then expect that the chaotic itinerant 
dynamics could be in two ways: chaotic transition dy- 
namics in time and in space. Not only the duration un- 
der the influence of one local attractor becomes chaotic, 
but so does the selection of the next-switching attractor. 
A detailed description of these systems including more 
than three local attractors will be reported elsewhere. 
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